rm(list = ls())

### R package
# # only need to install once
# devtools::install_github("cloudyr/pyMTurkR", dependencies=TRUE)
# devtools::install_github("Luwei-Ying/validateIt", dependencies=TRUE)
# load R package
library(validateIt)

### load the crowdsourced results
# word intrusion (WI)
wi_first <- vector() # first trial
wi_second <- vector() # second trial

load("topicrecords/WIstm10k1it.Rdata")
result <- read.csv("topicresults/WIstm10k1itresults1.csv", stringsAsFactors = FALSE)
wi_first[1] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])
result <- read.csv("topicresults/WIstm10k1itresults2.csv", stringsAsFactors = FALSE)
wi_second[1] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])

load("topicrecords/WIstm10k.Rdata")
result <- read.csv("topicresults/WIstm10kresults1.csv", stringsAsFactors = FALSE)
wi_first[2] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])
result <- read.csv("topicresults/WIstm10kresults2.csv", stringsAsFactors = FALSE)
wi_second[2] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])

load("topicrecords/WIstm50k.Rdata")
result <- read.csv("topicresults/WIstm50kresults1.csv", stringsAsFactors = FALSE)
wi_first[3] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])
result <- read.csv("topicresults/WIstm50kresults2.csv", stringsAsFactors = FALSE)
wi_second[3] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])

load("topicrecords/WIstm100k.Rdata")
result <- read.csv("topicresults/WIstm100kresults1.csv", stringsAsFactors = FALSE)
wi_first[4] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])
result <- read.csv("topicresults/WIstm100kresults2.csv", stringsAsFactors = FALSE)
wi_second[4] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])

load("topicrecords/WIstm500k.Rdata")
result <- read.csv("topicresults/WIstm500kresults1.csv", stringsAsFactors = FALSE)
wi_first[5] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])
result <- read.csv("topicresults/WIstm500kresults2.csv", stringsAsFactors = FALSE)
wi_second[5] <- as.numeric(evalResults(results = result, key = record, type = "WI")[[3]][1])

wi_pool <- (wi_first + wi_second)/2

# results so far should be as follows:
# wi_first <- c(0.212, 0.414, 0.512, 0.55, 0.518)
# wi_second <- c(0.252, 0.41, 0.476, 0.558, 0.47)
# wi_pool <- c(0.232, 0.412, 0.494, 0.554, 0.494)

# top 8 word set intrusion (T8WSI)
t8wsi_first <- vector() # first trial
t8wsi_second <- vector() # second trial

load("topicrecords/T8WSIstm10k1it.Rdata")
result <- read.csv("topicresults/T8WSIstm10k1itresults1.csv", stringsAsFactors = FALSE)
t8wsi_first[1] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])
result <- read.csv("topicresults/T8WSIstm10k1itresults2.csv", stringsAsFactors = FALSE)
t8wsi_second[1] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])

load("topicrecords/T8WSIstm10k.Rdata")
result <- read.csv("topicresults/T8WSIstm10kresults1.csv", stringsAsFactors = FALSE)
t8wsi_first[2] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])
result <- read.csv("topicresults/T8WSIstm10kresults2.csv", stringsAsFactors = FALSE)
t8wsi_second[2] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])

load("topicrecords/T8WSIstm50k.Rdata")
result <- read.csv("topicresults/T8WSIstm50kresults1.csv", stringsAsFactors = FALSE)
t8wsi_first[3] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])
result <- read.csv("topicresults/T8WSIstm50kresults2.csv", stringsAsFactors = FALSE)
t8wsi_second[3] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])

load("topicrecords/T8WSIstm100k.Rdata")
result <- read.csv("topicresults/T8WSIstm100kresults1.csv", stringsAsFactors = FALSE)
t8wsi_first[4] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])
result <- read.csv("topicresults/T8WSIstm100kresults2.csv", stringsAsFactors = FALSE)
t8wsi_second[4] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])

load("topicrecords/T8WSIstm500k.Rdata")
result <- read.csv("topicresults/T8WSIstm500kresults1.csv", stringsAsFactors = FALSE)
t8wsi_first[5] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])
result <- read.csv("topicresults/T8WSIstm500kresults2.csv", stringsAsFactors = FALSE)
t8wsi_second[5] <- as.numeric(evalResults(results = result, key = record, type = "T8WSI")[[3]][1])

t8wsi_pool <- (t8wsi_first + t8wsi_second)/2

# results so far should be as follows:
# t8wsi_first <- c(0.34, 0.526, 0.612, 0.654, 0.648)
# t8wsi_second <- c(0.34, 0.564, 0.62, 0.65, 0.612)
# t8wsi_pool <- c(0.340, 0.545, 0.616, 0.652, 0.630) 

# random 4 word set intrusion (R4WSI)
r4wsi_first <- vector() # first trial
r4wsi_second <- vector() # second trial

load("topicrecords/R4WSIstm10k1it.Rdata")
result <- read.csv("topicresults/R4WSIstm10k1itresults1.csv", stringsAsFactors = FALSE)
r4wsi_first[1] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])
result <- read.csv("topicresults/R4WSIstm10k1itresults2.csv", stringsAsFactors = FALSE)
r4wsi_second[1] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])

load("topicrecords/R4WSIstm10k.Rdata")
result <- read.csv("topicresults/R4WSIstm10kresults1.csv", stringsAsFactors = FALSE)
r4wsi_first[2] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])
result <- read.csv("topicresults/R4WSIstm10kresults2.csv", stringsAsFactors = FALSE)
r4wsi_second[2] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])

load("topicrecords/R4WSIstm50k.Rdata")
result <- read.csv("topicresults/R4WSIstm50kresults1.csv", stringsAsFactors = FALSE)
r4wsi_first[3] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])
result <- read.csv("topicresults/R4WSIstm50kresults2.csv", stringsAsFactors = FALSE)
r4wsi_second[3] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])

load("topicrecords/R4WSIstm100k.Rdata")
result <- read.csv("topicresults/R4WSIstm100kresults1.csv", stringsAsFactors = FALSE)
r4wsi_first[4] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])
result <- read.csv("topicresults/R4WSIstm100kresults2.csv", stringsAsFactors = FALSE)
r4wsi_second[4] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])

load("topicrecords/R4WSIstm500k.Rdata")
result <- read.csv("topicresults/R4WSIstm500kresults1.csv", stringsAsFactors = FALSE)
r4wsi_first[5] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])
result <- read.csv("topicresults/R4WSIstm500kresults2.csv", stringsAsFactors = FALSE)
r4wsi_second[5] <- as.numeric(evalResults(results = result, key = record, type = "R4WSI")[[3]][1])

r4wsi_pool <- (r4wsi_first + r4wsi_second)/2

# # results so far should be as follows:
# r4wsi_first <- c(0.258, 0.788, 0.84, 0.904, 0.71)
# r4wsi_second <- c(0.292, 0.806, 0.854, 0.92, 0.722)
# r4wsi_pool <- c(0.275, 0.797, 0.847, 0.912, 0.716)


### generate Figure 3
# define transparent colors for the figure
library("ggplot2")
forestgreen_t <- alpha("forestgreen", 0.35)
deepskyblue4_t <- alpha("deepskyblue4", 0.35)
darkorange2_t <- alpha("darkorange2", 0.35)
firebrick3_t <- alpha("firebrick3", 0.35)
mediumpurple3_t <- alpha("mediumpurple3", 0.35)
forestgreen <- "forestgreen"
deepskyblue4 <- "deepskyblue4"
darkorange2 <- "darkorange2"
firebrick3 <- "firebrick3"
mediumpurple3 <- "mediumpurple3"

# output plot in pdf
pdf("topicplots/figure3.pdf", width = 10, height = 7)
par(mgp = c(1.5, 0, 0), mar = c(2, 3, .7, .7))
plot(NULL,
     main = NA,
     ylim=c(0, 1.02), 
     xlim = c(0.6, 3.4), 
     ylab = "Proportion Correct",
     xlab = NA,
     cex.lab = 1.2,
     axes = F)
axis(side = 2, at = seq(0, 1, by = 0.2), col.ticks = NA, cex.axis = 1.2)
axis(side = 1, at = 1,
     labels = c('Word Intrusion'), 
     las = 1,
     col = NA,
     col.ticks = NA,
     cex.axis = 1.2)
axis(side = 1, at = 2,
     labels = c('Top 8\nWord Set Intrusion'),
     las = 1,
     col = NA,
     col.ticks = NA,
     cex.axis = 1.2)
axis(side = 1, at = 3,
     labels = c('Random 4\nWord Set Intrusion'),
     las = 1,
     col = NA,
     col.ticks = NA,
     cex.axis = 1.2)
legend(0.66, 1.06,
       c("Model 1: 10-topic, one iteration",
         "Model 2: 10-topic",
         "Model 3: 50-topic",
         "Model 4: 100-topic",
         "Model 5: 500-topic"),
       col = c(forestgreen, deepskyblue4, darkorange2, firebrick3, mediumpurple3),
       lty = 1,
       lwd = 3,
       cex = 1.2,
       bty = 'n')
# ----------------------------------------------------------------
segments(x0 = 0, y0 = 0.2, x1 = 1.5, y1 = 0.2, col = "gray", lty = 1, lwd = 2)
segments(x0 = 1.5, y0 = 0.25, x1 = 3.6, y1 = 0.25, col = "gray", lty = 1, lwd = 2)
# -------------------------------- WI ---------------------------------
# 10k1it
points(x = 0.7, y = wi_first[1], pch = 20, col = forestgreen_t, cex = 2)
segments(x0 = 0.7, y0 = wi_first[1]-1.96*sqrt(wi_first[1]*(1-wi_first[1])/500), x1 = 0.7, y1 = wi_first[1]+1.96*sqrt(wi_first[1]*(1-wi_first[1])/500), col = forestgreen_t, lwd = 3)
points(x = 0.72, y = wi_second[1], pch = 20, col = forestgreen_t, cex = 2)
segments(x0 = 0.72, y0 = wi_second[1]-1.96*sqrt(wi_second[1]*(1-wi_second[1])/500), x1 = 0.72, y1 = wi_second[1]+1.96*sqrt(wi_second[1]*(1-wi_second[1])/500), col = forestgreen_t, lwd = 3)
points(x = 0.75, y = wi_pool[1], pch = 20, col = forestgreen, cex = 2)
segments(x0 = 0.75, y0 = wi_pool[1]-1.96*sqrt(wi_pool[1]*(1-wi_pool[1])/1000), x1 = 0.75, y1 = wi_pool[1]+1.96*sqrt(wi_pool[1]*(1-wi_pool[1])/1000), col = forestgreen, lwd = 3)
# ----------------------------------------------------------------
# 10k
points(x = 0.85, y = wi_first[2], pch = 20, col = deepskyblue4_t, cex = 2)
segments(x0 = 0.85, y0 = wi_first[2]-1.96*sqrt(wi_first[2]*(1-wi_first[2])/500), x1 = 0.85, y1 = wi_first[2]+1.96*sqrt(wi_first[2]*(1-wi_first[2])/500), col = deepskyblue4_t, lwd = 3)
points(x = 0.87, y = wi_second[2], pch = 20, col = deepskyblue4_t, cex = 2)
segments(x0 = 0.87, y0 = wi_second[2]-1.96*sqrt(wi_second[2]*(1-wi_second[2])/500), x1 = 0.87, y1 = wi_second[2]+1.96*sqrt(wi_second[2]*(1-wi_second[2])/500), col = deepskyblue4_t, lwd = 3)
points(x = 0.9, y = wi_pool[2], pch = 20, col = deepskyblue4, cex = 2)
segments(x0 = 0.9, y0 = wi_pool[2]-1.96*sqrt(wi_pool[2]*(1-wi_pool[2])/1000), x1 = 0.9, y1 = wi_pool[2]+1.96*sqrt(wi_pool[2]*(1-wi_pool[2])/1000), col = deepskyblue4, lwd = 3)
# 50k
points(x = 1, y = wi_first[3], pch = 20, col = darkorange2_t, cex = 2)
segments(x0 = 1, y0 = wi_first[3]-1.96*sqrt(wi_first[3]*(1-wi_first[3])/500), x1 = 1, y1 = wi_first[3]+1.96*sqrt(wi_first[3]*(1-wi_first[3])/500), col = darkorange2_t, lwd = 3)
points(x = 1.02, y = wi_second[3], pch = 20, col = darkorange2_t, cex = 2)
segments(x0 = 1.02, y0 = wi_second[3]-1.96*sqrt(wi_second[3]*(1-wi_second[3])/500), x1 = 1.02, y1 = wi_second[3]+1.96*sqrt(wi_second[3]*(1-wi_second[3])/500), col = darkorange2_t, lwd = 3)
points(x = 1.05, y = wi_pool[3], pch = 20, col = darkorange2, cex = 2)
segments(x0 = 1.05, y0 = wi_pool[3]-1.96*sqrt(wi_pool[3]*(1-wi_pool[3])/1000), x1 = 1.05, y1 = wi_pool[3]+1.96*sqrt(wi_pool[3]*(1-wi_pool[3])/1000), col = darkorange2, lwd = 3)
# 100k
points(x = 1.15, y = wi_first[4], pch = 20, col = firebrick3_t, cex = 2)
segments(x0 = 1.15, y0 = wi_first[4]-1.96*sqrt(wi_first[4]*(1-wi_first[4])/500), x1 = 1.15, y1 = wi_first[4]+1.96*sqrt(wi_first[4]*(1-wi_first[4])/500), col = firebrick3_t, lwd = 3)
points(x = 1.17, y = wi_second[4], pch = 20, col = firebrick3_t, cex = 2)
segments(x0 = 1.17, y0 = wi_second[4]-1.96*sqrt(wi_second[4]*(1-wi_second[4])/500), x1 = 1.17, y1 = wi_second[4]+1.96*sqrt(wi_second[4]*(1-wi_second[4])/500), col = firebrick3_t, lwd = 3)
points(x = 1.2, y = wi_pool[4], pch = 20, col = firebrick3, cex = 2)
segments(x0 = 1.2, y0 = wi_pool[4]-1.96*sqrt(wi_pool[4]*(1-wi_pool[4])/1000), x1 = 1.2, y1 = wi_pool[4]+1.96*sqrt(wi_pool[4]*(1-wi_pool[4])/1000), col = firebrick3, lwd = 3)
# 500k 
points(x = 1.3, y = wi_first[5], pch = 20, col = mediumpurple3_t, cex = 2)
segments(x0 = 1.3, y0 = wi_first[5]-1.96*sqrt(wi_first[5]*(1-wi_first[5])/500), x1 = 1.3, y1 = wi_first[5]+1.96*sqrt(wi_first[5]*(1-wi_first[5])/500), col = mediumpurple3_t, lwd = 3)
points(x = 1.32, y = wi_second[5], pch = 20, col = mediumpurple3_t, cex = 2)
segments(x0 = 1.32, y0 = wi_second[5]-1.96*sqrt(wi_second[5]*(1-wi_second[5])/500), x1 = 1.32, y1 = wi_second[5]+1.96*sqrt(wi_second[5]*(1-wi_second[5])/500), col = mediumpurple3_t, lwd = 3)
points(x = 1.35, y = wi_pool[5], pch = 20, col = mediumpurple3, cex = 2)
segments(x0 = 1.35, y0 = wi_pool[5]-1.96*sqrt(wi_pool[5]*(1-wi_pool[5])/1000), x1 = 1.35, y1 = wi_pool[5]+1.96*sqrt(wi_pool[5]*(1-wi_pool[5])/1000), col = mediumpurple3, lwd = 3)
# ------------------------------- p-values -------------------------------
segments(x0 = 0.75, y0 = 0.1, x1 = 0.75, y1 = 0.19, col = "gray", lty = 3)
segments(x0 = 0.75, y0 = 0.1, x1 = 0.9, y1 = 0.1, col = "gray", lty = 3)
segments(x0 = 0.9, y0 = 0.1, x1 = 0.9, y1 = 0.38, col = "gray", lty = 3)
text(x = 0.83, y = 0.07, "0.000")
segments(x0 = 0.9, y0 = 0.6, x1 = 0.9, y1 = 0.46, col = "gray", lty = 3)
segments(x0 = 0.9, y0 = 0.6, x1 = 1.05, y1 = 0.6, col = "gray", lty = 3)
segments(x0 = 1.05, y0 = 0.6, x1 = 1.05, y1 = 0.54, col = "gray", lty = 3)
text(x = 0.98, y = 0.63, "0.000")
segments(x0 = 1.05, y0 = 0.4, x1 = 1.05, y1 = 0.45, col = "gray", lty = 3)
segments(x0 = 1.05, y0 = 0.4, x1 = 1.2, y1 = 0.4, col = "gray", lty = 3)
segments(x0 = 1.2, y0 = 0.4, x1 = 1.2, y1 = 0.51, col = "gray", lty = 3)
text(x = 1.13, y = 0.37, "0.008")
segments(x0 = 1.2, y0 = 0.66, x1 = 1.2, y1 = 0.6, col = "gray", lty = 3)
segments(x0 = 1.2, y0 = 0.66, x1 = 1.35, y1 = 0.66, col = "gray", lty = 3)
segments(x0 = 1.35, y0 = 0.66, x1 = 1.35, y1 = 0.54, col = "gray", lty = 3)
text(x = 1.28, y = 0.69, "0.000")
# -------------------------------- Task 2 ---------------------------------
# 10k1it
points(x = 1.7, y = t8wsi_first[1], pch = 20, col = forestgreen_t, cex = 2)
segments(x0 = 1.7, y0 = t8wsi_first[1]-1.96*sqrt(t8wsi_first[1]*(1-t8wsi_first[1])/456), x1 = 1.7, y1 = t8wsi_first[1]+1.96*sqrt(t8wsi_first[1]*(1-t8wsi_first[1])/456), col = forestgreen_t, lwd = 3)
points(x = 1.72, y = t8wsi_second[1], pch = 20, col = forestgreen_t, cex = 2)
segments(x0 = 1.72, y0 = t8wsi_second[1]-1.96*sqrt(t8wsi_second[1]*(1-t8wsi_second[1])/500), x1 = 1.72, y1 = t8wsi_second[1]+1.96*sqrt(t8wsi_second[1]*(1-t8wsi_second[1])/500), col = forestgreen_t, lwd = 3)
points(x = 1.75, y = t8wsi_pool[1], pch = 20, col = forestgreen, cex = 2)
segments(x0 = 1.75, y0 = t8wsi_pool[1]-1.96*sqrt(t8wsi_pool[1]*(1-t8wsi_pool[1])/1000), x1 = 1.75, y1 = t8wsi_pool[1]+1.96*sqrt(t8wsi_pool[1]*(1-t8wsi_pool[1])/1000), col = forestgreen, lwd = 3)
# 10k
points(x = 1.85, y = t8wsi_first[2], pch = 20, col = deepskyblue4_t, cex = 2)
segments(x0 = 1.85, y0 = t8wsi_first[2]-1.96*sqrt(t8wsi_first[2]*(1-t8wsi_first[2])/500), x1 = 1.85, y1 = t8wsi_first[2]+1.96*sqrt(t8wsi_first[2]*(1-t8wsi_first[2])/500), col = deepskyblue4_t, lwd = 3)
points(x = 1.87, y = t8wsi_second[2], pch = 20, col = deepskyblue4_t, cex = 2)
segments(x0 = 1.87, y0 = t8wsi_second[2]-1.96*sqrt(t8wsi_second[2]*(1-t8wsi_second[2])/500), x1 = 1.87, y1 = t8wsi_second[2]+1.96*sqrt(t8wsi_second[2]*(1-t8wsi_second[2])/500), col = deepskyblue4_t, lwd = 3)
points(x = 1.9, y = t8wsi_pool[2], pch = 20, col = deepskyblue4, cex = 2)
segments(x0 = 1.9, y0 = t8wsi_pool[2]-1.96*sqrt(t8wsi_pool[2]*(1-t8wsi_pool[2])/1000), x1 = 1.9, y1 = t8wsi_pool[2]+1.96*sqrt(t8wsi_pool[2]*(1-t8wsi_pool[2])/1000), col = deepskyblue4, lwd = 3)
# 50k
points(x = 2, y = t8wsi_first[3], pch = 20, col = darkorange2_t, cex = 2)
segments(x0 = 2, y0 = t8wsi_first[3]-1.96*sqrt(t8wsi_first[3]*(1-t8wsi_first[3])/500), x1 = 2, y1 = t8wsi_first[3]+1.96*sqrt(t8wsi_first[3]*(1-t8wsi_first[3])/500), col = darkorange2_t, lwd = 3)
points(x = 2.02, y = t8wsi_second[3], pch = 20, col = darkorange2_t, cex = 2)
segments(x0 = 2.02, y0 = t8wsi_second[3]-1.96*sqrt(t8wsi_second[3]*(1-t8wsi_second[3])/500), x1 = 2.02, y1 = t8wsi_second[3]+1.96*sqrt(t8wsi_second[3]*(1-t8wsi_second[3])/500), col = darkorange2_t, lwd = 3)
points(x = 2.05, y = t8wsi_pool[3], pch = 20, col = darkorange2, cex = 2)
segments(x0 = 2.05, y0 = t8wsi_pool[3]-1.96*sqrt(t8wsi_pool[3]*(1-t8wsi_pool[3])/1000), x1 = 2.05, y1 = t8wsi_pool[3]+1.96*sqrt(t8wsi_pool[3]*(1-t8wsi_pool[3])/1000), col = darkorange2, lwd = 3)
# 100k
points(x = 2.15, y = t8wsi_first[4], pch = 20, col = firebrick3_t, cex = 2)
segments(x0 = 2.15, y0 = t8wsi_first[4]-1.96*sqrt(t8wsi_first[4]*(1-t8wsi_first[4])/500), x1 = 2.15, y1 = t8wsi_first[4]+1.96*sqrt(t8wsi_first[4]*(1-t8wsi_first[4])/500), col = firebrick3_t, lwd = 3)
points(x = 2.17, y = t8wsi_second[4], pch = 20, col = firebrick3_t, cex = 2)
segments(x0 = 2.17, y0 = t8wsi_second[4]-1.96*sqrt(t8wsi_second[4]*(1-t8wsi_second[4])/500), x1 = 2.17, y1 = t8wsi_second[4]+1.96*sqrt(t8wsi_second[4]*(1-t8wsi_second[4])/500), col = firebrick3_t, lwd = 3)
points(x = 2.2, y = t8wsi_pool[4], pch = 20, col = firebrick3, cex = 2)
segments(x0 = 2.2, y0 = t8wsi_pool[4]-1.96*sqrt(t8wsi_pool[4]*(1-t8wsi_pool[4])/1000), x1 = 2.2, y1 = t8wsi_pool[4]+1.96*sqrt(t8wsi_pool[4]*(1-t8wsi_pool[4])/1000), col = firebrick3, lwd = 3)
# 500k
points(x = 2.3, y = t8wsi_first[5], pch = 20, col = mediumpurple3_t, cex = 2)
segments(x0 = 2.3, y0 = t8wsi_first[5]-1.96*sqrt(t8wsi_first[5]*(1-t8wsi_first[5])/500), x1 = 2.3, y1 = t8wsi_first[5]+1.96*sqrt(t8wsi_first[5]*(1-t8wsi_first[5])/500), col = mediumpurple3_t, lwd = 3)
points(x = 2.32, y = t8wsi_second[5], pch = 20, col = mediumpurple3_t, cex = 2)
segments(x0 = 2.32, y0 = t8wsi_second[5]-1.96*sqrt(t8wsi_second[5]*(1-t8wsi_second[5])/500), x1 = 2.32, y1 = t8wsi_second[5]+1.96*sqrt(t8wsi_second[5]*(1-t8wsi_second[5])/500), col = mediumpurple3_t, lwd = 3)
points(x = 2.35, y = t8wsi_pool[5], pch = 20, col = mediumpurple3, cex = 2)
segments(x0 = 2.35, y0 = t8wsi_pool[5]-1.96*sqrt(t8wsi_pool[5]*(1-t8wsi_pool[5])/1000), x1 = 2.35, y1 = t8wsi_pool[5]+1.96*sqrt(t8wsi_pool[5]*(1-t8wsi_pool[5])/1000), col = mediumpurple3, lwd = 3)
# ------------------------------- p-values -------------------------------
segments(x0 = 1.75, y0 = 0.2, x1 = 1.75, y1 = 0.3, col = "gray", lty = 3)
segments(x0 = 1.75, y0 = 0.2, x1 = 1.9, y1 = 0.2, col = "gray", lty = 3)
segments(x0 = 1.9, y0 = 0.2, x1 = 1.9, y1 = 0.5, col = "gray", lty = 3)
text(x =1.83, y = 0.17, "0.000")
segments(x0 = 1.9, y0 = 0.71, x1 = 1.9, y1 = 0.615, col = "gray", lty = 3)
segments(x0 = 1.9, y0 = 0.71, x1 = 2.05, y1 = 0.71, col = "gray", lty = 3)
segments(x0 = 2.05, y0 = 0.71, x1 = 2.05, y1 = 0.65, col = "gray", lty = 3)
text(x = 1.98, y = 0.74, "0.002")
# -------------------------------- R4WSI ---------------------------------
# 10k1it
points(x = 2.7, y = r4wsi_first[1], pch = 20, col = forestgreen_t, cex = 2)
segments(x0 = 2.7, y0 = r4wsi_first[1]-1.96*sqrt(r4wsi_first[1]*(1-r4wsi_first[1])/490), x1 = 2.7, y1 = r4wsi_first[1]+1.96*sqrt(r4wsi_first[1]*(1-r4wsi_first[1])/490), col = forestgreen_t, lwd = 3)
points(x = 2.72, y = r4wsi_second[1], pch = 20, col = forestgreen_t, cex = 2)
segments(x0 = 2.72, y0 = r4wsi_second[1]-1.96*sqrt(r4wsi_second[1]*(1-r4wsi_second[1])/500), x1 = 2.72, y1 = r4wsi_second[1]+1.96*sqrt(r4wsi_second[1]*(1-r4wsi_second[1])/500), col = forestgreen_t, lwd = 3)
points(x = 2.75, y = r4wsi_pool[1], pch = 20, col = forestgreen, cex = 2)
segments(x0 = 2.75, y0 = r4wsi_pool[1]-1.96*sqrt(r4wsi_pool[1]*(1-r4wsi_pool[1])/990), x1 = 2.75, y1 = r4wsi_pool[1]+1.96*sqrt(r4wsi_pool[1]*(1-r4wsi_pool[1])/990), col = forestgreen, lwd = 3)
# 10k  
points(x = 2.85, y = r4wsi_first[2], pch = 20, col = deepskyblue4_t, cex = 2)
segments(x0 = 2.85, y0 = r4wsi_first[2]-1.96*sqrt(r4wsi_first[2]*(1-r4wsi_first[2])/500), x1 = 2.85, y1 = r4wsi_first[2]+1.96*sqrt(r4wsi_first[2]*(1-r4wsi_first[2])/500), col = deepskyblue4_t, lwd = 3)
points(x = 2.87, y = r4wsi_second[2], pch = 20, col = deepskyblue4_t, cex = 2)
segments(x0 = 2.87, y0 = r4wsi_second[2]-1.96*sqrt(r4wsi_second[2]*(1-r4wsi_second[2])/500), x1 = 2.87, y1 = r4wsi_second[2]+1.96*sqrt(r4wsi_second[2]*(1-r4wsi_second[2])/500), col = deepskyblue4_t, lwd = 3)
points(x = 2.9, y = r4wsi_pool[2], pch = 20, col = deepskyblue4, cex = 2)
segments(x0 = 2.9, y0 = r4wsi_pool[2]-1.96*sqrt(r4wsi_pool[2]*(1-r4wsi_pool[2])/1000), x1 = 2.9, y1 = r4wsi_pool[2]+1.96*sqrt(r4wsi_pool[2]*(1-r4wsi_pool[2])/1000), col = deepskyblue4, lwd = 3)
# 50k
points(x = 3, y = r4wsi_first[3], pch = 20, col = darkorange2_t, cex = 2)
segments(x0 = 3, y0 = r4wsi_first[3]-1.96*sqrt(r4wsi_first[3]*(1-r4wsi_first[3])/500), x1 = 3, y1 = r4wsi_first[3]+1.96*sqrt(r4wsi_first[3]*(1-r4wsi_first[3])/500), col = darkorange2_t, lwd = 3)
points(x = 3.02, y = r4wsi_second[3], pch = 20, col = darkorange2_t, cex = 2)
segments(x0 = 3.02, y0 = r4wsi_second[3]-1.96*sqrt(r4wsi_second[3]*(1-r4wsi_second[3])/500), x1 = 3.02, y1 = r4wsi_second[3]+1.96*sqrt(r4wsi_second[3]*(1-r4wsi_second[3])/500), col = darkorange2_t, lwd = 3)
points(x = 3.05, y = r4wsi_pool[3], pch = 20, col = darkorange2, cex = 2)
segments(x0 = 3.05, y0 = r4wsi_pool[3]-1.96*sqrt(r4wsi_pool[3]*(1-r4wsi_pool[3])/1000), x1 = 3.05, y1 = r4wsi_pool[3]+1.96*sqrt(r4wsi_pool[3]*(1-r4wsi_pool[3])/1000), col = darkorange2, lwd = 3)
# 100k
points(x = 3.15, y = r4wsi_first[4], pch = 20, col = firebrick3_t, cex = 2)
segments(x0 = 3.15, y0 = r4wsi_first[4]-1.96*sqrt(r4wsi_first[4]*(1-r4wsi_first[4])/500), x1 = 3.15, y1 = r4wsi_first[4]+1.96*sqrt(r4wsi_first[4]*(1-r4wsi_first[4])/500), col = firebrick3_t, lwd = 3)
points(x = 3.17, y = r4wsi_second[4], pch = 20, col = firebrick3_t, cex = 2)
segments(x0 = 3.17, y0 = r4wsi_second[4]-1.96*sqrt(r4wsi_second[4]*(1-r4wsi_second[4])/500), x1 = 3.17, y1 = r4wsi_second[4]+1.96*sqrt(r4wsi_second[4]*(1-r4wsi_second[4])/500), col = firebrick3_t, lwd = 3)
points(x = 3.2, y = r4wsi_pool[4], pch = 20, col = firebrick3, cex = 2)
segments(x0 = 3.2, y0 = r4wsi_pool[4]-1.96*sqrt(r4wsi_pool[4]*(1-r4wsi_pool[4])/1000), x1 = 3.2, y1 = r4wsi_pool[4]+1.96*sqrt(r4wsi_pool[4]*(1-r4wsi_pool[4])/1000), col = firebrick3, lwd = 3)
# 500k 
points(x = 3.3, y = r4wsi_first[5], pch = 20, col = mediumpurple3_t, cex = 2)
segments(x0 = 3.3, y0 = r4wsi_first[5]-1.96*sqrt(r4wsi_first[5]*(1-r4wsi_first[5])/500), x1 = 3.3, y1 = r4wsi_first[5]+1.96*sqrt(r4wsi_first[5]*(1-r4wsi_first[5])/500), col = mediumpurple3_t, lwd = 3)
points(x = 3.32, y = r4wsi_second[5], pch = 20, col = mediumpurple3_t, cex = 2)
segments(x0 = 3.32, y0 = r4wsi_second[5]-1.96*sqrt(r4wsi_second[5]*(1-r4wsi_second[5])/500), x1 = 3.32, y1 = r4wsi_second[5]+1.96*sqrt(r4wsi_second[5]*(1-r4wsi_second[5])/500), col = mediumpurple3_t, lwd = 3)
points(x = 3.35, y = r4wsi_pool[5], pch = 20, col = mediumpurple3, cex = 2)
segments(x0 = 3.35, y0 = r4wsi_pool[5]-1.96*sqrt(r4wsi_pool[5]*(1-r4wsi_pool[5])/1000), x1 = 3.35, y1 = r4wsi_pool[5]+1.96*sqrt(r4wsi_pool[5]*(1-r4wsi_pool[5])/1000), col = mediumpurple3, lwd = 3)
# ------------------------------- p-values -------------------------------
segments(x0 = 2.75, y0 = 0.15, x1 = 2.75, y1 = 0.226, col = "gray", lty = 3)
segments(x0 = 2.75, y0 = 0.15, x1 = 2.9, y1 = 0.15, col = "gray", lty = 3)
segments(x0 = 2.9, y0 = 0.15, x1 = 2.9, y1 = 0.75, col = "gray", lty = 3)
text(x =2.83, y = 0.12, "0.000")
segments(x0 = 2.9, y0 = 0.92, x1 = 2.9, y1 = 0.84, col = "gray", lty = 3)
segments(x0 = 2.9, y0 = 0.92, x1 = 3.05, y1 = 0.92, col = "gray", lty = 3)
segments(x0 = 3.05, y0 = 0.92, x1 = 3.05, y1 = 0.88, col = "gray", lty = 3)
text(x = 2.98, y = 0.95, "0.000")
segments(x0 = 3.05, y0 = 0.7, x1 = 3.05, y1 = 0.8, col = "gray", lty = 3)
segments(x0 = 3.05, y0 = 0.7, x1 = 3.2, y1 = 0.7, col = "gray", lty = 3)
segments(x0 = 3.2, y0 = 0.7, x1 = 3.2, y1 = 0.88, col = "gray", lty = 3)
text(x = 3.13, y = 0.67, "0.008")
segments(x0 = 3.2, y0 = 0.97, x1 = 3.2, y1 = 0.94, col = "gray", lty = 3)
segments(x0 = 3.2, y0 = 0.97, x1 = 3.35, y1 = 0.97, col = "gray", lty = 3)
segments(x0 = 3.35, y0 = 0.97, x1 = 3.35, y1 = 0.77, col = "gray", lty = 3)
text(x = 3.28, y = 1, "0.000")
dev.off()